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We revisit the issue of cosmological parameter estimation in light of current and upcoming high- 
precision measurements of the cosmic microwave background power spectrum. Physical quantities 
which determine the power spectrum are reviewed, and their connection to familiar cosmological 
parameters is explicated. We present a set of physical parameters, analytic functions of the usual 
cosmological parameters, upon which the microwave background power spectrum depends linearly 
(or with some other simple dependence) over a wide range of parameter values. With such a set of pa- 
rameters, microwave background power spectra can be estimated with high accuracy and negligible 
computational effort, vastly increasing the efficiency of cosmological parameter error determination. 
The techniques presented here allow calculation of microwave background power spectra 10"" times 
faster than comparably accurate direct codes (after precomputing a handful of power spectra) . We 
discuss various issues of parameter estimation, including parameter degeneracies, numerical preci- 
sion, mapping between physical and cosmological parameters, and systematic errors, and illustrate 
these considerations with an idealized model of the MAP experiment. 

PACS numbers: 98.70.V, 98.80.C 



I. INTRODUCTION 

By January 2003, microwave maps of the full sky at 
0.2° resolution will be available to the world, the har- 
vest of the remarkable MAP satellite currently taking 
data 1^. The angular power spectrum of the tempera- 
ture fluctuations in these maps will be determined to high 
precision on angular scales from the resolution limit up to 
a dipole variation; this corresponds to about 800 statis- 
tically independent power spectrum measurements. Re- 
cent measurements have given a taste of the data to come, 
although covering much smaller patches of the sky and 
with potentially more serious systematic errors ||, 

The main science driving these spectacular technical 
feats is the determination of basic cosmological parame- 
ters describing our Universe, and resulting insights into 
fundamental physics; see |^, ^ for recent reviews and 
for a pedagogical introduction. The expected series 
of acoustic peaks in the microwave background power 
spectrum encode enough information to make possible 
the determination of numerous cosmological parameters 
simultaneously These parameters include the long- 
sought Hubble parameter Hq , the large-scale geometry of 
the Universe f2, the mean density of baryons in the Uni- 
verse rib, and the value of the mysterious but now widely 
accepted cosmological constant A, along with parameters 
describing the tiny primordial perturbations which grew 
into present structures, and the redshift at which the Uni- 
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verse reionized due to the formation of the first stars or 
other compact objects. While most of these parameters 
and other information will be determined to high preci- 
sion, one near-exact degeneracy and other approximate 
degeneracies exist between these parameters j^, |l^ . 

The exciting prospects of providing definitive answers 
to some of cosmology's oldest questions raise a poten- 
tially difficult technical issue, namely finding constraints 
on a large-dimensional parameter space. Given a set of 
data, what range of points in parameter space give mod- 
els with an acceptably good fit to the data? The answer 
requires evaluating a likelihood function at many points 
in parameter space; in particular, finding confidence re- 
gions in multidimensional parameter space requires look- 
ing around in the space. This is a straightforward pro- 
cess, in principle. But a parameter space with as many as 
ten dimensions requires evaluating a lot of models: a grid 
with a crude 10 values per parameter contains ten billion 
models. A direct calculation of this many models is pro- 
hibitive, even with very fast computers. While some of 
the parameters arc independent of the others (i.e. tensor 
modes), reducing the effective dimensionality, many of 
the parameters will be constrained quite tightly, requir- 
ing a finer sampling in that direction of parameter space. 
An additional problem with brute-force grid-based meth- 
ods is a lack of fiexibility: if additional parameters are 
required to correctly describe the Universe, vast amounts 
of recalculation must be done. Grid-based techniques 
have been used for analysis of cosmological parameters 
from microwave background 



[) but clearly 

this method is not fast or fiexible enough to deal with 
upcoming data sets adequately. 

A more sophisticated approach is to perform a search 
of parameter space in the region of interest. Relevant 
techniques are well known and have been applied to the 
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microwave background |T^, |T^, . Reliable estimates of 
the error region in cosmological parameter space can be 
obtained with random sets of around 10^ models, reduc- 
ing the computational burden by a factor of 1000 or more 
compared to the cruder grid search methods. On fast 
parallel computers it is possible to compute the power 
spectrum for a given model in a computation time on 
the order of a second jl^, making Monte Carlo error 
determinations feasible. For improved efficiency, useful 
implementations do not recalculate the entire spectrum 
at each point in parameter space, but rather use approx- 
imate power spectra based on smaller numbers of calcu- 
lated models 

Here we expand and refine this idea by presenting a set 
of parameters, functions of the usual cosmological param- 
eters, which reflect the underlying physical effects deter- 
mining the microwave background power spectrum and 
thus result in particularly simple and intuitive parame- 
ter dependences. Previous crude implementations of this 
idea applied to low- resolution measurements of the power 
spectrum ; the current and upcoming power spectrum 
measurements requires a far more refined implementa- 
tion. Our set of parameters can be used to construct 
computationally trivial but highly accurate approximate 
power spectra, and large Monte Carlo computations can 
then be performed with great efficiency. Additional ad- 
vantages of a physically-based parameter set are that the 
degeneracy structure of the parameter space can be seen 
much more clearly, and the Monte Carlo itself takes sig- 
nificantly fewer models to converge. 

Earlier Fisher-matrix approximations of parameter er- 
rors are a rough implementation of this general idea: if 
the the power spectrum varies exactly linearly with each 
parameter in some parameter set, and the measurement 
errors are gaussian random distributed and uncorrelated, 
then the likelihood function as a function of the parame- 
ters can be computed exactly from the partial derivatives 
of the model with respect to the parameters |8[ po| . 
Even if the linear parameter dependence does not hold 
throughout the parameter space, it will almost always be 
valid in some small enough region, being the lowest-order 
term in a Taylor expansion. If this region is at least as 
large as the resulting error region, then the Fisher matrix 
provides a self-consistent approximation for error deter- 
mination. The difficulty is in finding a suitable parameter 
set. Our aim is not necessarily to find a set of parame- 
ters which all have a perfectly linear effect on the power 
spectrum. Rather, more generally, we desire a set of pa- 
rameters for which the power spectrum can be approxi- 
mated very accurately with a minimum of computational 
effort. Then we can dispense with approximations of the 
likelihood function, as in the Fisher matrix approach, 
and directly implement a Monte Carlo parameter space 
search with a minimum of computational effort. This 
technique allows simple incorporation of prior probabil- 
ities determined from other sources of data, and if fast 
enough allows detailed exploration of potential system- 
atic biases in parameter values arising from numerous 



experimental and analysis issues such as treatment of 
foregrounds, noise correlations, mapmaking techniques, 
and model accuracy. 

We emphasize that with the accuracy of upcoming 
microwave background power spectrum measurements, 
the cosmological parameters will be determined precisely 
enough that their values will be significantly affected by 
systematic errors from assumptions made in the analysis 
pipeline and, potentially, numerical errors in evaluating 
theoretical power spectra. The measurements themselves 
may also be dominated by systematic errors. Extensive 
modelling of the impact of various systematic errors on 
cosmological parameter determination will be essential 
before any parameter determination can be considered 
reliable. This is the primary motivation for increasing 
the efficiency of parameter error analysis. We discuss 
this point in more detail in the last Section of the paper. 

The following Section reviews physical processes affect- 
ing the cosmic microwave background and defines our 
physical parameter set. The key parameter, which sets 
the angular scale of the acoustic oscillations, has been dis- 
cussed previously in similar contexts [|lO[ ^ ^ but has 
not been explicitly used as a parameter. The other pa- 
rameters describing the background cosmology can then 
be chosen to model other specific physical effects. Sec- 
tion [II displays how the power spectrum varies as each 



parameter changes with the others held fixed. The map- 
ping between our physical parameters and the conven- 
tional cosmological parameters immediately reveals the 
structure of degeneracies in the cosmological parameter 
space. Simple approximations for the effect of each phys- 
ical parameter are shown to be highly accurate for all 
but the largest scales in the power spectrum. We then 
determine the error region in parameter space for an ide- 
alized model of the MAP experiment in Sec. [V, com- 



paring with previous calculations, and test the accuracy 
of our approximate power spectrum calculation. Finally, 
Sec.0 discusses the potential speed of our method, like- 
lihood estimation, mapping between the sets of parame- 
ters, accuracy, and systematic errors. The first appendix 
summarizes an analysis pipeline going from a measured 
power spectrum to cosmological parameter error esti- 
mates; the second appendix details several numerical dif- 
ficulties with using the CMBFAST code for the calcula- 
tions in this paper, along with suggested fixes. 

Other recent attempts to speed up power spectrum 
computation through approximations include the dASH 
numerical package and interpolation schemes 
Our method has the advantage of conceptual simplicity 
and ease of use combined with great speed and high ac- 
curacy. 



II. COSMOLOGICAL PARAMETERS AND 
PHYSICAL QUANTITIES 

We consider the standard class of inflation-like cosmo- 
logical models, specified by five parameters determining 
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the background homogeneous spacetime (matter density 
f^mat, radiation density firad, vacuum energy density flj^, 
baryon density flf,, and Hubble parameter h), four pa- 
rameters determining the spectrum of primordial pertur- 
bations (scalar and tensor amplitudes As and At and 
power law indices n and tit), and a single parameter r de- 
scribing the total optical depth since reionization. We ne- 
glect additional complications such as massive neutrinos 
or a varying vacuum equation of state. This generic class 
of simple cosmological models appears to be a very good 
description of the Universe on large scales, and has the 
virtue of arising naturally in inflationary cosmology. Of 
course, any analysis of new data should first test whether 
this class of models actually provides an acceptably good 
fit to the data, and only then commence with determi- 
nations of cosmological parameters. 

What constraints do a measurement of the microwave 
background power spectrum place on this parameter 
space? A rough estimate assumes the likelihood of a par- 
ticular set of parameters is a quadratic function of the 
parameters (the usual Fisher matrix technique); better 
is a Monte-Carlo exploration of parameter space in the 
region of a best-fit model. For either technique to give 
reliable results, it is imperative to isolate the physical 
quantities which affect the microwave background spec- 
trum and understand their relation to the cosmological 
parameters. If a set of physical quantities which have 
essentially orthogonal effects on the power spectrum can 
be isolated, then extracting parameter space constraints 
becomes much more reliable and efficient. 

The five parameters describing the background cosmol- 
ogy induce complex dependences in the microwave back- 
ground power spectrum through multiple physical effects. 
Using the physical energy densities Slniat^i^, ^bh^, and 
flAh^, as has been done in numerous prior analyses, im- 
proves the situation but is still not ideal. The character- 
istic physical scale in the power spectrum is the angular 
scale of the first acoustic peak, so it is advantageous to 
use this scale as a parameter which can be varied inde- 
pendently of any other parameters. This angular scale is 
in turn determined by the ratio of the comoving sound 
horizon at last scattering (which determines the physi- 
cal wavelength of the acoustic waves) to the angular di- 
ameter distance to the surface of last scattering (which 
determines the apparent angular size of this yardstick). 
We first determine this quantity in terms of the cosmo- 
logical parameters, and then choose three other comple- 



mentary quantities. The analytic theory underlying any 
such choice of physical parameters has been worked out 
in detail (see ||, |2|, |2|l, and also H, |2|). 

The angular diameter distance to an object at a red- 
shift z is defined to be Da{z) = RpSkjr) / (1 + z) where, 
following the notation of Peacock , Rq is the (dimen- 
sionful) scale factor today and 

{sinr, f7 > 1; 
r! - 1; (1) 
sinhr, < 1; 
the Friedmann equation provides the connection between 
redshift z and coordinate distance r. With a dimension- 
less scale factor a and a Hubble parameter H = a/ a, the 
differential relation is 
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[(1 - rj)a2 + f^Aa^ + rimata + f^rad 
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Thus the relation between scale factor and coordinate 
distance becomes 



r = 

»i 

X 



1/2 



dx 
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(3) 



la [(1 - fl)x^ + ilAX^ + flmatX + f^rad 

and the angular diameter distance is just 

DA(a) = aH^'\n - l\-^^^Sk{r), (4) 

using the above expression for r. Here we have used the 
relation HqRq — \n ~ which follows directly from 

the Friedmann equation. 

The sound horizon is defined, in analogy to the particle 
horizon, as 



rsit) 



cjf) 
a{t') 



dt' 



(5) 



where Cs (t) is the sound speed of the baryon-photon fluid 
at time i; to a very good approximation, before decou- 
pling the sound speed is given by pj| 



= ^ (1 + 3pfc/4p^) \ 



(6) 
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Using this expression plus the equivalent expression to Eq. (y) connecting dt and da gives 

1 r dx 



rs{a) = 



(l + l^^x) ((1 - n)x^ + nAX'^ + anat2^ + ^ad) 



1/2- 



(7) 



The physical quantity relevant to the microwave back- ground power spectrum is 



^^^rr\ (8) 
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where the two functions are given by Eqs. and and 
a* is the scale factor at decoupUng. An accurate analytic 
fit for a, in terms of the cosmological parameters has 
been given by Hu and Sugiyama |2q]: 
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— - 1 = 1048 [l + 0.0012A(nbh^)-° '^^^] 

a* 

X [l+5l(f^mat/l')^^] , 

0.0783(f^6/i2)-°-238 [i + 39.5(f7fc/i2)0-763] 
0.560 [l + 21.1(f]b/i2)i-8i]^\ 



(9) 



This fit apphes to standard thermodynamic recombina- 
tion for a wide range of cosmological models. (Note that 
while the formula was constructed only for models with 
the standard value for flrad, the effect of a change of ra- 
diation density will come only through the redshift of 
matter-radiation equality, and thus ri,„at^^ can be re- 
placed by a factor proportional to ilmat/^^rad to account 
for a variation in the radiation density. The changes in 
are small anyway and have little impact on our analysis.) 
Since a, <C 1, the term proportional to J^a in Eg . (u h may 
be dropped and the integral performed, giving |25|| 
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where R = 3pb/4p-y is proportional to the scale factor. 

If two cosmological models are considered which both 
have adiabatic perturbations and the same value for A, 
to high accuracy their acoustic peaks will differ only in 
height, not position. This is an obvious advantage for 
constraining models. The peak positions, and thus A, 
will be extremely well constrained by power spectrum 
measurements. Other combinations which vary the peak 
heights only will be less well determined. 

Choosing four other parameters describing the back- 
ground cosmology must balance several considerations: 
(1) the new set of parameters must cover a sufficiently 
large region of parameter space; (2) the power spectrum 
should vary linearly or in some other simple way with the 
new parameters; (3) the new parameters should be nearly 
orthogonal in the cosmological parameter space; (4) the 
new parameters should correspond to the most impor- 
tant independent physical effects determining the power 
spectrum; (5) common theoretical prior constraints, like 
flatness or standard radiation, should be simple to im- 
plement by fixing a single parameter. No parameter set 
can satisfy all of these constraints perfectly. We use the 
following parameters, which provide a good balance be- 
tween these criteria: 



B 
V 

n 



^rad 



B is proportional to the baryon-photon density and thus 
determines the baryon driving effect on the acoustic os- 
cillations |2j, ^ . TZ is the matter-radiation density ratio 
at recombination, which determines the amount of early 
Integrated Sachs- Wolfe effect. V determines the late-time 
Integrated Sachs- Wolfe effect arising from a late vacuum- 
dominated phase, but otherwise represents a nearly exact 
degeneracy (sometimes called the "geometrical degener- 
acy" ) . M couples only to other small physical effects and 
is an approximate degeneracy direction. This choice of 
parameters is not unique, but this set largely satisfies the 
above criteria. 

Given values for A, B, V, TZ, and A4, they can be 
inverted to the corresponding cosmological parameters 
by rewriting the definition of A in terms of B, V,TZ, M, 
and h, then searching in h until the desired value for A 
is obtained, fif, and ^Ia then follow immediately, while 
ilinat and f^rad Can be obtained with a few iterations to 
determine a precise value for a* . A less efficient but more 
straightforward method we have implemented is simply 
to search the 5-dimensional cosmological parameter space 
for the correct values. 

The other cosmological parameters which affect the mi- 
crowave background power spectrum, aside from tensor 
perturbations, are the reionization redshift and the am- 
plitude and power law index of the scalar perturbation 
power spectrum. For reionization, we use the physical 
parameter 
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(12) 



the factor by which the microwave background 
anisotropics at small scales are damped due to Compton 
scattering by free electrons after the Universe is reion- 
ized. At large scales, the temperature fiuctuations are 
suppressed by a smaller amount. 

The primordial perturbation amplitude cannot be 
measured directly, but only the amplitude of the mi- 
crowave background fluctuations. Some care must be 
taken in the definition of the amplitude of the microwave 
fiuctuations so that the other physical parameters are 
not significantly degenerate with a simple change in nor- 
malization. For example, one common normalization is 
some weighted average of the C;'s over the smallest I 
values, say 2 </< 20 111. This is known as "COBE- 
normalization" because this is the range of scales probed 
by COBE, which tightly constrains the normalization of 
the temperature fluctuations at these scales . But al- 
though such a normalization is useful for interpreting the 
COBE results, it is a bad normalization to adopt when 
probing a much larger range in /-space. The reason is 
that the lowest multipoles have a significant contribu- 
tion from the Integrated Sachs- Wolfe effect: if the total 
matter density changes significantly between two models, 
the resulting low I multipole moments will also change. 
If the normalization is fixed at low I, then the result is 
that the two models will be offset at high / values. As a 
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concrete example, consider two models which differ only 
in V, with A, B, M. and TZ held fixed. The only phys- 
ical difference between these two models is a difference 
in the structure growth rate at low redshifts, and in par- 
ticular the fluctuations at the last scattering surface are 
identical. But if these two models are COBE normal- 
ized, then every C/ for Z = 20 is offset between the two 
models. Defining the normalization as an average band 
power over a wider range of / values is better but not 
ideal. We advocate defining a normalization parameter 
S by the amplitude of the perturbations on the scale of 
the horizon at last scattering, corrected by the amount of 
small-scale suppression in the C/ 's due to reionization. In 
practice, for models with the same values of n, this nor- 
malization essentially corresponds to the amplitude of the 
acoustic oscillations at I values higher than the first sev- 
eral peaks. Note that increasing Z while holding S and 
the other physical parameters fixed keeps the microwave 
background power spectrum the same at high I while en- 
hancing the signal at low I. The first few peaks cannot be 
used for normalization because they have been subjected 
to driving by the gravitational potential as they cross the 
horizon p7[ , so their amplitudes vary significantly with 
B and TZ, as shown in the following Section. 

The primordial power spectrum of density perturba- 
tions is generally taken to be a power law, P{k) oc fc", 
with n — 1 corresponding to the scale-invariant Harrison- 
Zeldovich spectrum. Making the approximation that k 
and I have a direct correspondence, the effect of n on the 
microwave background power spectrum can be modelled 
as 



(13) 



which is a good approximation for power law power spec- 
tra. A departure from a power law, characterized by a 
p2| , can be represented in a similar way. Note that since 
the dependence is exponential, a linear extrapolation is 
never a good approximation over the entire interesting 
range 2 < I < 3000. The choice of Iq is arbitrary: chang- 
ing Iq simply gives a different overall normalization. 

Finally, note that the amplitude of the tensor mode 
power spectrum should be used as a separate parameter, 
not the common choice of the tensor-scalar ratio. The 
tensor modes contribute at comparatively large angular 
scales; they are significant only for I < 100 given cur- 
rent rough limits on the tensor amplitude. Varying the 
tensor amplitude independently automatically leaves the 
normalization parameter S fixed. 

In the rest of the paper, we refer to the parameters A, 
B, V, TZ, M, S, and Z as "physical parameters", as op- 
posed to the usual "cosmological parameters" Jlf,, finiat, 
firad, f^vac, ^, quadrupolc-based normalization), 

and Zr- 



III. POWER SPECTRA 

We now present the dependence of the temperature 
fluctuation power spectrum on these physical parame- 
ters. As a fiducial model, we choose a standard cosmol- 
ogy with n = 0.99, A = 0.7, r^mat = 0.29, fit ^ 0.04, 
h — 0.7, n — 1, and full reionization at redshift Zr — 7, 
with standard neutrinos. We compute power spectra for 
the fiducial model and for evaluating the various numeri- 
cal derivatives using Seljak and Zaldarriaga's CMBFAST 
code 1 33, 34, 3^. Note that the model is slightly open 
rather than flat to compensate for numerical instabili- 
ties in the CMBFAST code; see the discussion in the 
Appendix. The corresponding physical parameters are 
Ao = 0.0106, Bo = 0.0196, Vo = 0.338, Mo = 0.154, and 
TZo = 3.252. Since the dependence of the spectrum on 
S is trivial, its definition is arbitrary; it can be taken, 
for example, as the amplitude of a particular high peak, 
or as the equivalent scalar perturbation amplitude As at 
the scale which the power law n is defined for a model 
with no reionization. We do not consider tensor pertur- 
bations in this paper, which can be analyzed separately 
by the same techniques; the tensor perturbations are sim- 
pler because they depend on fewer parameters than the 
scalar perturbations. None of the following results de- 
pend qualitatively on the fiducial model chosen. 

The top panel of Fig. |l] displays power spectra as the 
parameter A varies while the other physical parameters 
are held fixed. It is clear that the effect of A is only to 
determine the overall angular scale, except for the small 
Integrated Sachs- Wolfe effect for I < 200. For the bot- 
tom panel, the ?-axis of each power spectrum is rescaled 
to give Ci as a function of lA/Ao, and then the percent 
difference between the rescaled model and the fiducial 
model is plotted. The small variation of the scaled power 
spectrum between around I = 50 and I — 200 is due to 
the early ISW effect and can be accurately modelled as a 
linear variation, while between I = 2 and I = 50 the vari- 
ation with A is more complicated due to the late-time 
ISW contribution. For I > 200, the scaled power spectra 
match to within 1%; the residual discrepancy is likely 
due to numerical inaccuracy. Since A will be tightly con- 
strained by the peak positions, this numerical error will 
not affect the determination of A, but may contribute a 
systematic error to the determination of other parame- 
ters at the 1% level. 

Figure |^ displays the power spectrum variation with 
B, keeping the other parameters fixed. This plot clearly 
shows the well-known baryon signature of alternating 
peak enhancement and suppression. The variation of a 
particular multipole Ci with B is highly linear for I > 50 
for variations of up to 20%; for all multipoles, the depen- 
dence is linear for variations of up to 10%. A figure of 
merit for a linear approximation to the power spectrum 
is whether the errors in the approximation are negligible 
compared to the cosmic statistical error at a given multi- 
pole. The bottom panel of Fig. ^ shows just how good the 
linear approximation is for several selected I values. The 
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FIG. 1: Temperature power spectra for the fiducial cosmological model (top panel, heavy line), plus models varying the 
parameter A upwards by 10% (curve shifts to lower I values) and downwards by 10% (curve shifts to higher I values) while 
keeping B, M, V , TZ, S and n fixed. The bottom panel displays the fractional error between the fiducial model and the other 
two models with i-axes rescaled by the factor A/Ai. 



horizontal axis plots the fractional change in the parame- 
ter B. The vertical axis plots the ratio of the change in C; , 
Ci{B) — Ci{Bo), to the statistical error from the cosmic 
variance at that multipole, ai = Ci{Bq)^/{2/{21 + 1)). 
(Any measurement of the multipole Ci is subject to a 
statistical error at least as large as this cosmic variance 
error, due to the finite number of modes sampled on the 



sky.) The linear dependence remains valid even near the 
"pivot" points between a peak and a trough. Thus a lin- 
ear extrapolation using computed numerical derivatives 
is highly accurate at reproducing the B dependence; the 
error in such an approximation appears to be dominated 
by systematic numerical errors in computing the power 
spectrum, and funny behavior of some of the multipoles 
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FIG. 2: Temperature power spectra for the fiducial cosmological model (top panel, heavy line), plus models varying the 
parameter B upwards by 25% (higher first peak) and downwards by 25% (lower first peak) while keeping A, A4, V, TZ, S, and 
n fixed. The bottom panel displays C; as a function of B, for / = 50, 100, 200, 500, and 1000. The horizontal axis shows the 
fractional change in B, while the vertical axis gives the change in Ci as a fraction of the cosmic variance at that multipole. 



for small variations in B is surely due to systematic er- 
rors. 

Figure ^ shows the variation of the power spectrum 
with respect to the parameter TZ, keeping the others 
fixed, with the corresponding dependence of specific mul- 
tipoles. For smaller values of TZ, matter-radiation equal- 
ity occurs later, so the Universe is less accurately de- 



scribed as matter-dominated at the time of last scatter- 
ing, leading to an increased amplitude of the first few 
peaks from the Integrated Sachs- Wolfe effect at early 
times. The dependence on TZ is not quite as linear as 
for B, but still the spectrum varies highly linearly with 
TZ for parameter variations of 10% for / > 30. Figure ^ 
shows the same for the parameter A4 , which is an approx- 
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FIG. 3: Same as Fig. except varying the parameter TZ while keeping the others fixed. Larger TZ increases the peak heights. 
The substantial glitch in the bottom panel is due to systematic errors in CMBFAST. 



imate degeneracy direction: the power spectrum changes 
noticably only in the neighborhood of the first peak. If 
the standard three massless neutrino species is imposed 
as a prior, A4 is fixed. 

Variation of V holding the others fixed is nearly an ex- 
act physical degeneracy j|, 0, broken only by the late- 
time ISW effect at low multipole moments displayed in 
Fig. ^ This variation at small / values is not well approx- 
imated by simple linear extrapolation; accurate approx- 



imation must rely on more sophisticated schemes, but 
the total contribution to constraining cosmological mod- 
els will have comparatively little weight due to the large 
cosmic variance. At higher / values, the power spectra 
should be precisely degenerate, and the extent to which 
they are not is a measure of numerical inaccuracies. Note 
that the inaccuracies are in the form of systematic offsets 
at levels smaller than a percent. When analyzing data, 
the variation in V can simply be set to zero for all I val- 
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FIG. 4: Same as Fig. ^ except varying the parameter M while keeping the others fixed. Larger M increases the first peak 
height. Again, the glitch in the bottom panel is due to systematic errors in CMBFAST. 



ues higher than the first acoustic peak. If a flat universe 
is imposed as a prior condition, the parameter V is then 
fixed. 

Aside from initial conditions, reionization is the other 
physical effect which can have a significant impact on the 
power spectrum. Note that the effective normalization 
we use depends on Z, so that as Z varies, the power spec- 
trum retains the same amplitude at high /. As discussed 
in the previous section, for large I the effect of chang- 



ing the total optical depth is effectively just a change in 
the overall amplitude. Thus for our set of parameters, 
variations in Z while holding the other parameters fixed 
change the power spectrum only at the largest scales. 
For this reason, we do not consider reionization further 
here. Note that for polarization, reionization will gener- 
ate a new peak in the polarization power spectrum at low 
I, and this effect may require more complicated analytic 
modelling than a simple linear extrapolation ||3^ . 
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FIG. 5: Varying the the parameter V, keeping the others fixed, for low I (top) and higher I (center) values; note the scales of 
the vertical axes. Only the lowest multipoles have any significant variation, arising from the Integrated Sachs- Wolfe effect at 
late times; the bottom panel shows the dependence of the / = 5, 10, and 20 multipoles on V. A linear approximation is rough 
but reasonable; higher-order approximations will model this dependence better. At higher I, the variation between the curves 
is a measure of the numerical accuracy of the code generating the Ci curves. 
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FIG. 6: A comparison of computed (solid) and approximated (dashed) power spectra for different values of tlie scalar index n, 
for n = 0.8 and n = 1.2 (dashed), where the approximated values have been calulated by applying Eq. ( p^ ) to the n — 1 model 
(dotted). Note the displayed variation in n is roughly three times larger than the constraint MAP will impose (see Fig. Bt). 



Finally, we show the effect of changing the spectral in- 
dex n of the primordial power spectrum in Fig. ^. Equa- 
tion ([l3| ) gives a very good approximation to the effect 
of changing n; the largest errors are at small I and in 
the immediate region of the first few peaks and troughs, 
amounting to less than 10% of the cosmic variance for 
any individual multipole. 

We have now displayed simple numerical or analytic 
approximations for the effect of all parameters with a 
sizeable impact on the cosmic microwave background 
power spectrum, for I > 50 and for parameter variations 
which are in the range to which the parameters will be re- 
stricted by upcoming microwave background maps. For 
a model which differs from a fiducial model in several 
different parameters, it is in principle ambiguous in what 
order to do the various extrapolations and approxima- 
tions. For example, should a reionization correction be 
applied before or after the A parameter extrapolation? In 
practice, however, most of the parameter ranges involved 
will be quite small, particularly for the A parameter, and 
the order in which all of the approximations are applied 
to the fiducial spectrum is essentially irrelevant. We test 
the validity of the various approximations for arbitrary 



models below. 

We also note that while we have demonstrated explic- 
itly the validity of linear extrapolations and other simple 
approximations in computing deviations from a particu- 
lar fiducial power spectrum, we have not shown explic- 
itly that this is true for any fiducial cosmological model. 
But current measurements point towards a Universe well- 
described by a model close to the fiducial one considered 
here, and in the unlikely event that the actual cosmology 
is significantly different, the various approximations are 
simple to check explicitly for any underlying model. 

Some of the linear extrapolations lose accuracy for the 
smallest / values. It is not clear whether any simple ap- 
proximation is sufficient for reproducing the power spec- 
trum for small I values, because a variety of different 
effects contribute in varying proportions. While these 
large scales have significant cosmic variance and do not 
have very much statistical weight compared to the rest 
of the spectrum, several parameters (V, Z, T, ny) have 
significant effects only at small Z, and accurate limits on 
these parameters will require approximating these mo- 
ments. One possible approach is to use quadratic or 
higher-order extrapolations. We have not explored the 
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efficacy of such approximations, and they may require a 
more accurate numerical code, but promise to be a good 
solution. A more complicated approach would be to de- 
velop analytic fitting formulas. Steps in this direction 
have been taken by Durrer and collaborators |^, but 
their solutions rapidly lose accuracy for I > 10, and a va- 
riety of physical effects are neglected which contribute at 
the few percent level or more. Another possibility is an 
expanded numerical approximation, where the effects of 
our parameters are fit as linear or quadratic, but instead 
of fitting the total C; dependence, the various separate 
terms (i.e., quadratic combinations of Sachs- Wolfe, Inte- 
grated Sachs-Wolfe, doppler, and acoustic effects) con- 
tributing to each C/ are modelled independently. In the 
estimates below, we will not consider tensor perturba- 
tions and fix the parameter Z, while V is a degenerate 
direction; we thus sidestep the issue of approximating the 
power spectrum for low I, and do not discuss the issue 
further here. 



IV. ERROR REGIONS FOR COSMOLOGICAL 
PARAMETERS 

Given some data, we want to determine the error re- 
gion in parameter space corresponding to some certain 
confidence level in fitting the data. This must be done by 
looking around in parameter space in the vicinity of the 
best-fitting model. With a fast computation of the theo- 
retical power spectrum for different points in parameter 
space in hand, a straightforward Metropolis algorithm 
can be used to construct a Monte-Carlo exploration of 
the parameter space. Refined techniques with the label 
of Markov Chains have recently been applied to the mi- 
crowave background [ p^ |3^ . Our approximate power 
spectrum evaluation renders these techniques highly effi- 
cient. 

As a demonstration, we consider estimated power spec- 
trum measurement errors drawn from the X21+1 proba- 
bility distribution derived in Ref. We assume that 
the Universe is actually described by a particular fiducial 
cosmological model. With the simple assumption that 
the measurement errors for each C/ are uncorrelated, the 
likelihood of a cosmological model being consistent with 
a measurement of the fiducial model becomes trivial to 
compute. Then we map the likelihood around the fidu- 
cial model via a Markov Chain of points in parameter 
space, starting with the fiducial model and moving to a 
succession of random points in the space using a sim- 
ple Metropolis-Hastings algorithm. Sophisticated con- 
vergence tests can be applied to the chain to determine 
when the distribution of models in the Monte Carlo has 
converged to a representation of the underlying likelihood 
function. The chain of models should be constructed in 
the physical parameter space: the roughly orthogonal ef- 
fects of the physical parameters on the power spectrum 
leads to an efficient Monte Carlo sampling of the space. 
For each model in the chain, the cosmological parame- 



ters of the model are computed and stored along with 
the physical parameters. 

The commonly-used Fisher matrix approximation |^ 
expands the likelihood function in a Taylor series around 
the most likely model and retains only the quadratic 
term. With Gaussian uncorrelated statistical errors, the 
Fisher matrix approximation becomes exact whenever 
the parameter dependence of the Ci 's is exactly linear in 
all parameters. Since the power spectrum is nearly linear 
in most of our physical parameter set, we can check our 
Monte Carlo results by comparing with the Fisher matrix 
likelihood computed with the physical variables. 

As a simple illustration, we determine the error region 
corresponding to a simplified model of the power spec- 
trum which will be obtained by the MAP satellite, cur- 
rently collecting data. MAP's highest frequency channel 
has a gaussian beam with a fuU-width-half-max size of 
about 0.21 degrees. MAP will produce a full-sky map of 
the microwave background at this resolution with a sen- 
sitivity of around 35 /iK per 0.3 square-degree pixel. We 
neglect eventual sky cuts and assume all Ci estimates 
are uncorrelated. We consider models with only scalar 
perturbations and fixed reionization. The physical pa- 
rameter space is thus A, B, A4, TZ, V and n. The other 
parameters we have fixed (tensor perturbations, Z) af- 
fect the power spectrum only at small I, and will have 
little effect on the physical parameters considered. We as- 
sume the underlying cosmology is described by the fidu- 
cial model of the previous section. Fundamental physical 
constraints have been enforced: fib < fimat and all den- 
sities must be positive. Additionally, we only consider 
models with fivac < 2 and discard any models which 
are not invertible from the physical to the cosmological 
parameters (these amount to a handful of models with 
extreme parameter values). 

Figure |^ shows the anticipated MAP error contours in 
the physical parameter space, extracted from a Markov 
Chain of 3 x lO'' models. The full 7-dimensional like- 
lihood region is projected onto all pairs of parameters. 
The parameter V is essentially unconstrained by MAP, 
while A4 has a 1-a error of around 30%. Measurements 
extended to smaller scales will not significantly improve 
constraints on these two parameters: they represent true 
physical degeneracies in the model space. The other pa- 
rameters are well constrained. Approximate 1-a errors 
are: 0.5% for A, 5% for 7^, 3% for B, 2% for S, and 3% 
for n. The error regions are elliptical, because the like- 
lihood region in these parameters is well-approximated 
by a quadratic form. The residual correlations involving 
S and n will be lifted by measurements out to higher I 
values; the parameters are largely uncorrelated in their 
power spectrum effects. Figure |8| shows a Fisher matrix 
estimate of the same likelihood. The two sets of likeli- 
hood contours are very close, as expected since we have 
shown explicitly that the power spectrum varies nearly 
linearly in the physical variables over a parameter re- 
gion consistent with MAP-quality data. The most no- 
table exception is the variable A, which does not give 
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FIG. 7: Parameter error contours anticipated for the MAP satellite in the physical parameter space. In each panel, the entire 
likelihood region has been projected onto a particular plane. The axis refer to the fractional variation for each parameter with 
respect to the fiducial one, e.g. 5A = {A — ^fld)/^fld. The contours show 68%, 95%, and 99% confidence levels. 



a linear variation, but it is so well constrained that it is 
effectively fixed when considering its impact on determin- 
ing the other parameters. The contours will also differ 
slightly due to the nonlinear behavior of n, but the gen- 
eral excellent agreement demonstrates that our Monte 
Carlo technique works as expected. 

In Fig. H, we replot Fig. |^ in the cosmological parame- 
ter space. This procedure is trivial, since we calculate the 



cosmological parameter set for each model. The resulting 
contours are deformed, reflecting the non-linear relation 
between the physical and cosmological parameter sets. 
Clearly, a Fisher matrix approximation in the cosmolog- 
ical variables, as has been done in numerous analyses, 
will give significantly incorrect error regions ||l0|| . Note, 
however, that for MAP-quality data, the Fisher matrix 
errors in the physical variables projected into the cosmo- 
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FIG. 8: The same as Fig. except using the Fisher matrix approximation to the actual likelihood. The contours in the two 
cases are very similar. 



logical variables gives a highly accurate result, Fig. |lO|, in 
contrast to the claims in Ref . . The Fisher matrix ap- 
proximation in the physical parameter space can be used 
to further increase the efhciency of error region evalua- 
tion, which may be useful in the cases where the compu- 
ation of the error regions is dominated by evaluation of 
the likelihood for a given model from actual data. If the 
likelihood evaluation is at least as efficient as the power 
spectrum evaluation, then the increased efficiency from a 



Fisher matrix approximation may not result in dramatic 
increases in computational speed, since in both cases the 
physical parameter space must be sampled and parame- 
ter conversion performed at each point to construct the 
likelihood contours in the space of the cosmological pa- 
rameters. 

It is useful to plot error regions in both sets of param- 
eters. Physical parameters represent the actual physical 
quantities which are being directly constrained by the 
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FIG. 9: The likelihood distribution in Fig. ^ mapped into the cosmological parameter space. 



measurements, and the error regions are nearly elliptical. 
The cosmological parameters are useful theoretically and 
for comparisons with other data sources, i.e. large-scale 
structure or supernova standard candles. One advantage 
of determining the error region via a Monte Carlo tech- 
nique is that prior constraints on cosmological parame- 
ters can be implemented simply. 

When the Fisher matrix approximation is computed in 
the physical parameter space, it is highly accurate. This 
does not actually help much, however: if information 



about cosmological parameters is desired, transforming 
from the physical to the cosmological parameter space 
requires some kind of sampling of the likelihood region. 
The most efficient way to do this is via Monte Carlo, so 
the additional computational overhead in performing a 
Monte Carlo likelihood evaluation of the physical vari- 
able likelihood instead of a Fisher matrix approximation 
is negligible. 

To test the validity of our numerical power spectrum 
approximations, we choose 1000 models at random from 
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FIG. 10: The likelihood distribution in Fig. ^ mapped into the cosmological parameter space. These contours, derived from a 
Fisher matrix approximation, are very close to those in Fig. S. 



the above Monte Carlo set of cosmological models. We 
compute each numerically using CMBFAST, and then 
compare the computed and approximated models. The 
two are compared using the statistic 



X 



(14) 



where {(jf^^^)"^ is the variance of C; in our model MAP 
data. For MAP, the measured power spectrum will 



consist of approximately 800 independent multipole mo- 
ments, so an approximate power spectrum which differs 
from the actual model power spectrum by the MAP er- 
ror at each multipole will have a value for the statistic 
of around 800. Figure p] shows the error distribution 
for our subset of models. The statistic for each model 
is plotted against the "distance" in physical parameter 
space of the model from the fiducial model, where the 
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distance measure is the Fisher matrix: 

d=[(P-Po)F(P-Po)]'/' 



1 r 



(15) 



where P is a vector of physical parameters correspond- 
ing to the model considered, Pq is the vector of fiducial- 
model parameters, and 



E 



dCi dCi 



(16) 



is the Fisher matrix. This is roughly equivalent to mea- 
suring the distance along each axis in the physical pa- 
rameter space in units of the l-cr error interval in that 
parameter direction. The solid-circle points are models 
with 2 < < 4, for which the worst-approximated 
models disappear. The difference has two possible ori- 
gins: either linear extrapolation is not valid over a large 
enough range in the variables TZ and M to handle these 
models, or CMBFAST does not provide accurate compu- 
tations in this range. The second possibility seems more 
likely, since the code contains an explicit trap prevent- 
ing computations for these unusual numbers of neutrino 
species. 

For essentially all other models, the error averaged over 
all multipoles is small compared to the MAP statistical 
errors, and the error distribution is roughly independent 
of the distance in parameter space, showing that our nu- 
merical approximations are valid. Furthermore, most of 
the total error tends to come from a small number of mul- 
tipoles, indicating that systematic errors in CMBFAST 
dominate the total difference between our approxima- 
tions and the direct numerical calculations. Given this 
fact, we forego a more thorough characterization of the 
errors in our approximate spectrum computations until a 
more accurate code is available, but the approximations 
are likely accurate enough for reliable parameter analy- 
sis. It is unclear at this point whether the errors in our 
approximations or the systematic errors in CMBFAST 
will have a greater impact on the inferred error region. 
(Several other public CMB Boltzmann codes are in fact 
available, but none are as general as CMBFAST, allow- 
ing for arbitrary curvature, vacuum energy, and neutrino 
species, all of which are required for the analysis pre- 
sented here.) 



V. DISCUSSION AND CONCLUSIONS 

Error regions similar to Fig. ^ have been previously 
constructed for model microwave background experi- 
ments [^0[ The remarkable point about our calcu- 
lation is that the entire error region, constructed from 
power spectra for 3 x lO'' points in parameter space, has 
been computed in a few seconds of time on a desktop 
computer. In fact, since our power spectra can be com- 
puted with a few arithmetic operations per multipole mo- 
ment, the calculation time might not even be dominated 
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FIG. 11; A plot of the difference between approximate and nu- 
merical power spectra versus distance from the fiducial model 
in parameter space, for 1000 cosmological models drawn ran- 
domly from the distribution in Fig. fzl The solid circles are 
for models with 2 < A^^ < 4. 



by computing model power spectra, but rather by com- 
puting the likelihood function or by converting between 
the cosmological and physical parameters. 

For actual non-diagonal covariance matrices describing 
real experiments, the computation of the likelihood will 
dominate the total computation time. Gupta and Heav- 
ens have recently formulated a method for computing an 
approximate likelihood with great efficiency |Q, based 
on finding uncorrelatcd linear combinations of the power 
spectrum estimates |4^. Each parameter corresponds 
to a unique combination, so the likelihood is calculated 
with a handful of operations. Moreover, the method is 
optimized so the parameter estimation can be done with 
virtually no loss of accuracy. Using this technique, like- 
lihood estimates for realistic covariance matrices will be 
roughly as efficient as our power spectrum estimates. 

The conversion from the physical parameters to cos- 
mological ones requires evaluation of several numerical 
integrals. The integrands will be smooth functions, and 
this step can likely be made nearly as efficient as the 
power spectrum evaluation, although we have not imple- 
mented an optimized routine for parameter conversion. 
We make the following rough timing estimate: each mul- 
tipole of the power spectrum can be approximated by a 
handful of floating point operations, and an approximate 
likelihood and the parameter conversion both will plausi- 
bly require similar computational work. So we anticipate 
on the order of 10 floating point operations per multipole 
for each of 1000 multipoles. On a 1 gigaflop machine, this 
results in evaluation of around 10^ models per second, so 
a Markov Chain with tens of thousands of models can be 
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computed in under a second. Of course, such a compu- 
tation is trivially parallelizable, and can easily be sped 
up by a factor of tens on a medium-size parallel machine. 
In contrast, other state-of-the-art parameter estimation 
techniques |l^, , dominated by the power spectrum 
calculation, compute roughly one model per second, re- 
quiring several hours for a lO'* point computation. 

This great increase in computational speed is highly 
useful. All upcoming microwave background measure- 
ments and resulting parameter estimates will be domi- 
nated by systematic errors, both from measurement er- 
rors when observing the sky and processing errors in the 
data analysis pipeline. The only way to discern the ef- 
fect of these errors is through modelling them, requiring 
a determination of parameters numerous times. The mi- 
crowave background has the potential to constrain pa- 
rameters at the few percent level. But a variation of this 
size in any given parameter will result in changes in Ci 's 
of a few percent, which is significantly smaller than the 
errors with which each Ci will be measured. This means 
that small systematic effects distributed over many CPs 
can bias derived parameter values by amounts larger than 
the formal statistical errors on the parameter. Exhaus- 
tive simulation of a wide range of potential systematics 
will be required before the accuracy of highly precise pa- 
rameter determinations can be believed, and our tech- 
niques make such investigations far faster and easier. 

Along with increased speed, our approximation meth- 
ods also promise high accuracy. A careful error evalua- 
tion is hindered by systematic errors in CMBFAST; our 
numerical approximations give better accuracy than the 
predominant numerical code for parameter ranges gen- 
erally larger than the region which will be allowed by 
the MAP data. The physical set of variables presented 
here clarify the degeneracy structure of the power spec- 
trum, clearly displaying degenerate and near-degenerate 
directions in the parameter space. For example, with the 
physical parameters presented here, it is clear that any 
tensor mode contribution, which affects only low-Z mul- 
tipoles, will have negligible impact on the error contours 
in Fig. 0. 

Bond and Efstathiou used a principle component anal- 
ysis to claim that uncertainties in the tensor mode con- 
tributions would be the dominant source of error for all 
of the cosmological parmeters, and further claimed that 
as a result, the Fisher matrix approximation will over- 
estimate the errors in other apparently well-determined 
parameters such as our B |jl^. These results are clearly 
valid only for data which is significantly less constraining 
than the MAP data will be. Essentially, their approxi- 
mate degeneracy trades off variations in TZ and which 
change the heights of the first few peaks, with variations 
in Us, which also impact the first few peaks, while holding 
A fixed. From our analysis, it is clear that this degen- 
eracy only holds approximately, since TZ and B have a 
significant impact only on the first three peaks, while Us 
affects all of the peak heights. Indeed, a careful exami- 
nation of Fig. 3 in shows that the power spectra are 



only approximately degenerate, with peak height differ- 
ences which are significant for MAP. The addition of ten- 
sor modes allows a rough match between the two power 
spectra at low I values, but for high / values the discrep- 
ancy in peak heights is already evident at the third peak 
and will become larger for higher peaks, since the values 
of Us in the two models are so different. 

Gravitational lensing makes a negligible difference in 
the linearity of the parameter dependences. We have not 
included polarization in this analysis, but the three other 
polarized CMB power spectra can be handled in the same 
way as the temperature case using standard techniques 

We have shown that a proper choice of physical param- 
eters enables the microwave background power spectrum 
to be simply approximated with high accuracy over a sig- 
nificant region of parameter space. This region will be 
large enough for analyses of data from the MAP satellite, 
although larger regions could be stitched together using 
multiple reference models. Our approximation scheme re- 
quires first computing numerical derivatives of the power 
spectrum multipoles with respect to the various parame- 
ters, but this must be done only once and then computing 
further spectra is extremely fast (in the neighborhood of 
10^ models per second or more on common computers). 
The error determination techniques in this paper require 
the numerical evaluation of only tens of power spectra 
rather than thousands or millions, so speed of the power 
spectrum code will not be of paramount concern. On 
the other hand, we emphasize that for any high-precision 
parameter estimate, even small systematic errors in com- 
puting C;'s can lead to biased parameter estimates com- 
parable to the size of the statistical errors. These two 
considerations argue for a significant revision of CMB- 
FAST (CMBSLOW, perhaps?) or construction of other 
independent codes which focus on overall accuracy and 
stability of derivatives with respect to parameters rather 
than on computational speed. Such a code, combined 
with the estimation techniques in this paper and efficient 
likelihood evaluation methods will provide a highly 
efficient and reliable way to constrain the space of fun- 
damental cosmological parameters. 
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APPENDIX A: A MODEL PIPELINE FOR 
PARAMETERS ANALYSIS 

In the interests of completeness, in this Appendix we 
outHne an analysis pipeline for cosmological parameters 
and their errors. Our basic data set will be a high- 
resolution microwave background map, from which has 
been extracted an angular power spectrum (i.e., the C;'s) 
and a covariance matrix for the C/'s. Given this data 
set, plus potentially relevant data from other sources, we 
need to find the best-fit cosmological model and the error 
contours in the parameter space. 

In the case of an ideal, full-sky map with only sta- 
tistical errors and uniform sky coverage, the covariance 
matrix reduces to a multiple of the identity matrix, but 
in general, for most experiments the coviance matrix may 
contain significant off-diagonal terms. Even for the MAP 
satellite, which will come close to a diagonal covariance 
matrix, the need to mask out the galactic plane induces 
some correlations between nearby multipoles. The first 
step in the analysis is to construct a likelihood calculator: 
what is the likelihood that the measured data represents 
some particular theoretical power spectrum? This is in 
general a computationally hard problem, since it involves 
inverting an N x N matrix, where N is the number of 
multipoles in the power spectrum. However, the prob- 
lem becomes easier if the covariance matrix is close to 
diagonal. Heavens and Gupta have recently formulated 
a method for computing an approximate likelihood with 
great efficiency [Q, based on an eigenvector decompo- 
sition of the covariance matrix They have demon- 
strated a likelihood calculator which takes an input the- 
oretical power spectrum and outputs a likelihood based 
on simulated measurements and covariances in much less 
than a second of computation time. Further efficiency 
improvements are desirable, as this is likely to be the 
piece of the analysis pipeline which dominates the total 
time. 

The likelihood calculator can also incorporate prior 
probabilities obtained from other data sources, or from 
physical considerations like A > 0. Since none of the fol- 
lowing steps depend on the form of the likelihood func- 
tion, the inclusion of other kinds of data can be handled 
in a completely modular way. 

The second element needed is a pair of routines to 
convert between the physical parameter space (.4, B, 
M, V) and the cosmological parameter space {^mat, ^b, 
Ha, h). To go from the cosmological parameters to the 
physical parameters simply involves evaluating the inte- 
gral needed in Eq. ^ The reverse direction can be done 
straightforwardly by rewriting A in terms of S, A4, V, 
TZ, and h and then numerically finding the value of h 
which gives the needed value of A. The total variation 
of A with h is a smooth function, and the root-finding 
step can be done with a minimal number of integral eval- 
uations. This procedure does not guarantee a solution, 
however, in extreme regions of parameter space. A more 
general (but slower) method is to search the cosmolog- 



ical parameter space for the model which best fits the 
given set of physical parameters, and then project this 
model back into the physical parameter space to check 
the quality of the fit. 

Once a likelihood calculator and parameter conversion 
routines are in hand, we advocate the following proce- 
dure for parameter determination. First, choose a fidu- 
cial model which gives a good fit to the measured power 
spectrum. This can largely be done by eye, using the pa- 
rameter dependences displayed in the Figures in Sec. IH. 
Note that the fiducial model does not need to be the for- 
mal best-fit model, but only needs to be close enough to 
the best-fit model so that models throughout the error 
region will be accurately approximated by the extrapo- 
lations and approximations in Sec. HI. 



Next, check the validity of the fiducial model. Is the fit 
of the model with the data adequate? It is certainly pos- 
sible that the simple space of cosmological models con- 
sidered here may not contain the actual universe. For 
example, some admixture of isocurvature initial condi- 
tions might be important, or the initial power spectrum 
might not be well described by a power law. If the best-fit 
model is not a good fit to the data, it makes no sense to 
proceed with any further parameter analysis! The model 
space must be enlarged or modified to include models 
which provide a good fit to the data. 

If the fiducial model is adequate, then compute nu- 
merical partial derivatives of each C; with respect to B, 
M, and TZ. Also, write numerical routines to compute 
the power spectrum approximately for Z < 50. Even 
though CMBFAST or similar codes run very efficiently 
when computing only 2 < I < 50, this time will still 
greatly dominate the total computation time for param- 
eter error determination, so an efficient analytic estimate 
is essential. Note that an alternate possibility is to sim- 
ply ignore the lowest multipoles, which will be a small 
portion of the data for a high-resolution map, at the ex- 
pense of losing any leverage on the parameters V, Z, and 
the tensor power spectrum. Still, the other parameters 
will be determined much more precisely, so it may be 
possible to fit the rest of the parameters using I > 50 
data only, and then constrain just the V-Z-tensor pa- 
rameter space from the I < 50 data independently. Now 
it is possible to estimate the C/ 's highly efficiently for any 
model within a sizeable parameter space region around 
the fiducial model. 

Now perform a Markov-chain Monte Carlo jl^ in the 
physical parameter space. Each point in parameter space 
requires evaluation of one approximate power spectrum, 
one likelihood, and one conversion to cosmological pa- 
rameters. The Monte Garlo will converge efficiently since 
the physical parameters are nearly orthogonal in their 
effects on the power spectrum. Sophisticated conver- 
gence diagnostics are available for determining when the 
Markov chain of models has sufficiently sampled the error 
region 

Finally, once the chain of models in parameter space 
has been computed, extract the best-fit model and error 
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contours from the distribution of models. By displaying 
two-dimensional likelihood plots for all parameter combi- 
nations, a reliable picture of the shape of the likelihood 
region in the whole multidimensional parameter space 
can be visualized. 

Upcoming data will be good enough that the statistical 
errors, given by these likelihood contours, will be smaller 
than the systematic errors, at least in some parameter 
space directions. Sources of systematic error include in- 
strumental effects, modelling of foreground emission, ap- 
proximation of the covariance matrix, inaccuracies in the 
theoretical models, and inaccuracies in the approxima- 
tions in this paper. Before we can confidently claim to 
have determined the cosmological parameters and their 
errors, these systematics must be modelled and investi- 
gated. Using the techniques outlined here, where a com- 
plete Monte Carlo can be accomplished in seconds, in- 
vestigation of systematics is easily practical, instead of a 
computational challenge. 



APPENDIX B: SOME NUMERICAL 
CONSIDERATIONS WITH THE CMBFAST 
CODE 

Scljak and Zaldarriaga's CMBFAST code ||, H has 
become the standard tool for computing the microwave 
background temperature and polarization power spectra 
for the inflation-type cosmological models discussed in 
this paper. Hundreds of researchers have employed it 
for tasks ranging from predictions of highly speculative 
cosmological models to analysis of every microwave back- 
ground measurement of the past few years. Its compre- 
hensive treatment of a variety of physical effects and its 
public availability have made it the gold standard for mi- 
crowave background analysis. 

The original aim of the CMBFAST code was fast eval- 
uation of the power spectra. Previous codes had simply 
evolved the complete Boltzmann hierarchy of equations 
up to the maximum desired value of /. This resulted in a 
coupled set of around 3/ linear differential equations, and 
the /-values of interest might be 1000 or greater. The task 
of evolving thousands of coupled equations with oscilla- 
tory solutions led to codes which required many hours 
to run on fast computers. The CMBFAST algorithm 
evolves the source terms first, which involve angular mo- 
ments only up to / = 4; the same Boltzmann hierarchy 
can be used, but with a cutoff of / ~ 10. The rest of 
the C;'s can then be obtained by integrating this source 
term against various Bessel functions. The Bessel func- 
tions can be precomputed or approximated |^6| indepen- 
dently for each I. Great time savings result since (1) it is 
not necessary to compute the power spectrum for every Z, 
but rather only, e.g., every 50th I with the rest obtained 
from interpolation; (2) the time steps for the Boltzmann 
integration are not controlled by the need to resolve all 
of the oscillations in the Bessel functions, and thus the 
Boltzmann integrator can be significantly more efficient. 



It also may be possible to employ approximations to the 
integrals of the source times the Bessel functions with an 
additional gain in speed, although CMBFAST does not 
implement this. The code offers a speed improvement 
of a factor of 50 over Boltzmann hierarchy codes for flat 
universes, although the advantage is smaller for models 
which are not flat or have a cosmological constant. 

CMBFAST is nominally accurate at around the one 
percent level, although no systematic error analysis has 
ever been published. While the numerical error in any 
given multipole Ci is relatively small, the derivatives of 
Ci with respect to the various cosmological parameters 
are significantly less accurate. At the time the code was 
written, no microwave background experiments were of 
sufficient precision to warrant worrying about inaccura- 
cies at this level, and the great utility of the code has been 
its combination of speed and accuracy. But we have now 
entered the era of high-precision measurements which will 
place significant constraints on many parameters simul- 
taneously, and as this paper has demonstrated, efficiently 
prosecuting this analysis requires a code of high precision 
and stable behavior with respect to parameter variations. 
We have found that the unmodified version of CMBFAST 
"out of the box" is not sufficiently accurate for the calcu- 
lations presented here, but these shortcomings are par- 
tially compensated by the following procedures. 

First, CMBFAST uses slightly different algorithms for 
flat and non-flat background spacetime geometries. In 
the flat case, the necessary spherical Bessel functions 
are precomputed and called from a file; in the nonflat 
case, the hyperspherical Bessel functions are computed 
via an integral representation as they are needed. As 
a result, the line-of-sight integrals over the product of 
source term and Bessel function have different partitions 
in the two cases, and the numerical values differ at the 
percent level. Also, the code has a flag which forces 
evaluation of the flat case whenever |f2 — 1| < 0.001. If 
the flducial model is flat, then as parameters are varied 
which change the geometry, many individual multipoles 
Ci exhibit discontinuities at the point in parameter space 
where the code switches from the flat to nonflat evalua- 
tion scheme. These discontinuities can result in signifl- 
cant errors when evaluating Ci derivatives with respect 
to parameters. Also, even if the derivative is correctly 
evaluated, the offset between the two cases can result in 
a slightly biased error estimate. 

A careful flx of this problem would involve insuring 
that the integral partitions used in the two cases are de- 
termined consistently. Alternately, the flat-space code 
could be rewritten to numerically integrate the neces- 
sary Bessel functions as in the non-flat case, since the 
two methods do not have a signiflcant difference in com- 
putational speed. A simpler practical solution which we 
have employed in this paper is to force the code always to 
use the non-flat integration routine. This can be accom- 
plished by narrowing the tolerance at which the code uses 
the non-flat integration (say to 1| < 10~^), combined 
with always using a non-flat flducial model (shifting ilmat 
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FIG. 12: Derivatives of Ci with respect to B, calculated by 
CMBFAST with AZ = 10, from 2%, 5%, and 10% finite dif- 
ferences in B. The smoothest curve is for the 10% variation. 

by 10~^ will never result in a statistically significant dif- 
ference when fitting a given data set, due to cosmic vari- 
ance). 

Second, in varying the B parameter, adjacent acoustic 
peaks move in the opposite direction. Between the peaks, 
therefore, is a pivot point I for which Ci remains fixed as 
B is varied, and the nearby C;'s will vary only slightly. 
As CMBFAST is written, it computes C/ at I values sep- 
arated by 50; the complete spectrum is then obtained by 
splining between the calculated values. Splines are rather 
stiff, and small changes in regions where the power spec- 
trum has significant curvature (i.e. around the acoustic 
peaks and troughs) can noticeably change the overall fit 
to the power spectrum in the regions between. In partic- 
ular, the position of the "pivot" ^-value in the spectrum 
sometimes shifts a bit, which gives a spurious irregular 
dependence on B of the few C;'s in this region. 

This can be addressed by a simple modification to 
CMBFAST forcing the explicit evaluation of Ci for more 
I values. We have used a 51 of 10, which has eliminated 



the spurious behavior. Of course, this results in a signif- 
icant increase in the computational time needed for an 
individual model. A more sophisticated approach would 
be to use a comparable number of / values as the original 
code, but distribute them preferentially in the /-ranges 
corresponding to the peaks and troughs where the cur- 
vature of the power spectrum is largest. This is easily 
implemented using the variables defined in this paper, 
since the peak positions are completely determined by A 
alone. Using a finer grid in I also reproduces the peak 
positions with higher accuracy. 

Third, as discussed in the body of the paper, the power 
spectra must be normalized according to a physical cri- 
terion which does not depend on non-linear behavior at 
small I values. This is accomplished by simply compiling 
CMBFAST with the UNNORM option. In this case, no 
normalization to COBE is done; instead, the amplitude 
of the potential 5* is taken to be unity at the scale of the 
horizon at the time of last scattering. This normalization 
is fixed by the primordial amplitude of scalar perturba- 
tions, often denoted as As- We define a physical variable 
S corresponding to the microwave background amplitude 
by also including the effect of differing amounts of reion- 
ization so that the power spectrum at large I remains 
fixed. 

Some residual numerical problems remain with CMB- 
FAST. We have found, for example, that in the region 
around / — 350 (between the first and second acoustic 
peaks for our fiducial model), the power spectrum ex- 
hibits spurious oscillatory behavior with an amplitude of 
a few percent, if computed directly for every I instead of 
every 10th or 50th / (see Fig. |2|). This likely arises from 
either the k or rj integral over an oscillatory integrand not 
being evaluated on a fine enough grid. As a result, di- 
rectly calculated C/ values in this region of the parameter 
space effectively have non- negligible random errors. This 
problem is normally hidden by the spline smoothing with 
every 50th I computed. When the number of C; evalu- 
ations is increased to every 10th I, the power spectrum 
shows a distinct glitch in this region, and the parameter 
dependence is completely unreliable here. The spline fits 
are forced into more oscillations to include the computed 
points. More accurate numerical integrations will prob- 
ably solve this problem, at the expense of a significant 
speed reduction. 
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